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The thermodynamic and dynamical behavior of a gas of hard 
disks in a narrow channel is studied theoretically and numerically. 
Using a virial expansion we find that the pressure and collision fre- 
quency curves exhibit a singularity at a channel width corresponding 
to twice the disk diameter. As expected, the maximum Lyapunov 
exponent is also found to display a similar behavior. At high density 
these curves are dominated by solid-like configurations which are dif- 
ferent from the bulk ones, due to the channel boundary conditions. 

I. INTRODUCTION 

The thermodynamic behavior of hard-disk systems has been a subject of con- 
siderable interest in recent years [1]. In particular, the freezing transition taking 
place in these systems has been analysed by computer simulations for cases where 
the aspect ratio of the simulation box is close to one. Furthermore, the dynamical 
stability of such systems has been studied in detail and Lyapunov spectra have 
been computed [2,3]. The maximum Lyapunov exponent was found to exhibit 
a maximum in the density regime characteristic for the fluid-to-solid transition 
[4]. In a restricted geometry, such as a narrow channel, one expects other in- 
teresting features to emerge when the width of the channel becomes comparable 
to the disk's size. Studying narrow channels also enables one to examine long- 
wavelength perturbed states associated with the small positive Lyapunpov expo- 
nents found in the tangent-space dynamics of hard-disk and hard-sphere systems 
[2,3]. This provided the motivation for us to re- investigate the thermodynamic 
and dynamical properties of a hard-disk gas in a narrow channel. 
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It is found that for narrow channels the pressure exhibits a singularity at 
a channel width equal to twice the diameter of the disks. Narrower channels 
do not allow particles to pass each other. At the same width a singularity is 
also observed for the collision frequency. In addition, the maximum Lyapunov 
exponent and the Kolmogorov- Sinai entropy display a non-trivial dependence on 
the channel width in this regime. 

The model we study consists of N hard disks in a two-dimensional box with 
side lengths L^.Ly and aspect ratio A = Ly/Lj.. We vary the width of the box, 
keeping the volume V = L^Ly and the particle density p = N/V constant. Very 
small aspect ratios A are considered, for which the box resembles a narrow 
channel. For most of our work we use periodic boundary conditions in both direc- 
tions. Exact thermodynamic properties are known for extremely narrow channels 
in this case Ly < [5,6]. Recently, the transport properties of two particles in 
a square periodic box, A = 1, have been studied by computer simulation [7]. The 
case of two particles with reflecting boundary conditions in both directions has 
also received attention, both by molecular dynamics [8] and analytical approaches 
[9]. It was found that for certain aspect ratios the compressibility becomes neg- 
ative. We find analogous results also for partly-refiecting boundaries, which are 
refiecting for the walls parallel to the long direction of the box, x, and periodic 
for the walls parallel to y. 

Throughout, reduced units are used for which the disk diameter a, the particle 
mass m and the kinetic energy per particle, K/N, are unity. There is no potential 
energy in this case, and the total energy E is identical to the kinetic energy K. 
As a consequence, the temperature T is just an irrelevant parameter and is fixed 
to unity. Also Boltzmann's constant k is taken unity. 

The paper is organized as follows. In Section II we consider the low-density 
limit. In Section II A we present an analytical derivation of the pressure and 
show that the singularity taking place at L^, = 2 already exists in this limit. In 
addition, we obtain an analytical expression for the collision frequency which 
shows a behavior very similar to that of the pressure. Results of numerical 
simulations of these quantities and of the Lyapunov spectra are presented in 
Section II B. A numerical study of the high-density case is presented in Section 
III. Full Lyapunov spectra for periodic systems with periodic boundaries are 
discussed in Section IV. We close with a brief discussion of the results in Section 
V. 
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II. LOW-DENSITY LIMIT FOR SYSTEMS WITH PERIODIC 

BOUNDARIES 



A. Virial approach to the equation of state and collision frequency 

The pressure in the low-density hmit is obtained by the leading terms in 
the virial expansion as applied to the narrow-channel geometry. Let q{Ly) be 
the excluded volume of a single disk in the channel. The phase-space volume 
associated with placing N disks in the channel is given by 

n{N) ^V{V - q) . . .{V - {N - l)q), (1) 

which, to leading order in the density, is reduced to 

Here, v — V/N is the specific volume. For low density the specific entropy of the 
gas is given by 

where k is Boltzmann's constant. 

For Ly > 2 the excluded area is simply given by the area of a disk of radius 1, 
namely q{Ly > 2) — tt. To leading order in the density, the pressure P is given 

by 

Pv ds ^ TT 

and is thus independent of the channel width Ly. This is the expected leading- 
order expression of the virial expansion of a gas of hard disks. 

For 1 < Ly < 2, however, the excluded area is reduced due to the interference 
with the disk image resulting from the periodic boundary condition in the y 
direction (see Fig. 1). Simple geometrical considerations yield 

q{Ly) ^26 + sm{2e), (5) 

where 

Ly = 2sm{9). (6) 

The excluded volume q{Ly) depends on Ly. Taking the derivative of the entropy 
with respect to the volume and keeping the aspect ratio A = Ly/L^ constant, 
one finds to leading order in the density 
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FIG. 1. The excluded area associated with a disk in a channel of width Ly < 2 for 
the case of periodic boundary conditions. The area q{Ly) is given by that of the union 
of the disk and its translationally-displaced image within the channel. 



Note that the pressure curve exhibits a square root singularity as one approaches 
Ly = 2 from below. 

The singularity at Ly = 2 does not signify a phase transition in the usual sense 
as it is not a collective phenomenon. In fact, the singularity exists even for finite 

and it explicitly shows up in the second term of the virial expansion. This is 
different from the usual liquid-solid transition in the bulk, which is obtained only 
after summing over all orders of that expansion. The singularity in the narrow 
channel is a consequence of the fact that the available volume for the disks is a 
singular function of the control parameter Ly at Ly = 2. Higher order terms in 
the virial expansion are expected to exhibit a singularity at other integer channel 
widths, Ly = 3,4... . Thus singularities in the pressure are expected to show up 
at these widths at higher densities. 

We now turn to the evaluation of the collision frequency z/2. The virial ex- 
pression for the pressure, appropriately modified for two-dimensional hard disks, 
is given by [10] 



where the sum is over all collisions taking place during the time interval r. Here, 
the relative position vector of two colliding particles i and j at the time of the 




(7) 




(8) 



c 
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collision, Tc = — r^-, satisfies Ir^ = ex, and Vc = — Vj is the corresponding 
relative pre-coUision velocity. This may be written as 

PV 1 

- 1 - T^UF^^Ap (9) 



NkT 2kT 

where Ap is the average momentum transfer in a collision. For a large aspect 
ratio A the collision frequency is given by [11] 

U2 = 2n'/\kTf\a/v)g{a), (10) 

where g{r) is the pair distribution function. At contact, r = a, and in the low- 
density hmit one has g{a) = 1. Combining this result with the Eqs. (4) and (9), 
one finds for the average momentum transfer 

Ap = (v^/2)\/fcT. (11) 

This expression may also be derived directly from the Maxwell-Boltzmann veloc- 
ity distribution. It is expected to hold also for narrow channels. The reason is 
that in the narrow channel the velocity distribution is still isotropic and is given 
by the Maxwell-Boltzmann distribution (see Section 11 B). Therefore, one can 
use Eq. (11), together with Eq. (9), to obtain the collision frequency for narrow 
channels: 

Vkf V^\kT J ^ ' 



This expression is vahd both below and above Ly — 2. 

It is expected that in such a highly anisotropic system the pressure tensor is 
not isotropic. In the low-density limit one can evaluate the diagonal components 
of the pressure tensor, P^x and Pyy, using again the expression for the entropy in 
Eq. (3). One finds 



P — T 

^ XX ■'-'x 



Pyy ~ ^y 



ds 



dLj^^ v-{l/2)q{Lyy 
ds\ _v-{l/2)Ly^/4^^ 



dLyJr v-{l/2)q 



(13) 



Furthermore, one has 



fds\ 1 r ds \ I , f ds 



\dvj^ 2 ^\dLj^^ 2 '\dLy 

and, hence, P = {P^x + Pyy)/'^- 

It is evident that in the low-density limit the pressure tensor is isotropic for 
Ly > 2. It becomes anisotropic for Ly < 2. At higher densities the two pressure- 
tensor components differ even above Ly — 2 (see Section III). The square-root 
singularity of the total pressure P at Ly = 2 originates from Pyy. The other 
component, Pxx, exhibits a weaker singularity. 
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B. Numerical simulations 



Here we carry out numerical simulations for the computaton of the pressure, 
the collision frequency and the Lyapunov exponents. The pressure is evaluated 
from the impulsive version of the virial theorem Eq. (8) [12,13]. The collision 
frequency of a single particle, 1^2, is obtained from the total number of collisions 
per unit time, divided by N/2. For the Lyapunov spectrum we use the algorithm 
outlined in Refs. [4] and [2]. For the periodic-houndarj case the linear momentum 
is conserved both in x and y directions, and, as expected, altogether 6 Lyapunov 
exponents vanish due to the conservation of energy, momentum, center of mass 
(only in tangent space), and the regularity of the dynamics in the phase-flow 
direction. For the partly-reflecting boundary case only four of the exponents 
vanish due to a lack of momentum and center-of-mass conservation in the y 
direction. The Kolmogorov- Sinai (KS) entropy, hxs, which measures the rate of 
information gain on the initial conditions due to the time-reversible dynamics, 
is given by the sum of the positive exponents [14]. Since its dependence on the 
channel width is qualitatively similar to that of the maximum Lyapunov exponent 
Ai wc concentrate on the latter in the following. 

In Figure 2 we compare the theoretical curves for the pressure and colli- 
sion frequency, Eqs. (4), (7) and (12), with the numerical results. In this low- 
density regime an excellent agreement is obtained. Also shown in this figures is 
the channel-width dependence of the maximum Lyapunov exponent. It shows a 
shoulder at the critical width Ly = 2. The tangent-space perturbation associated 
with Ai has been shown to be localized in space [3]. This exponent is closely 
connected to the collision frequency and, as a consequence, has a similar qual- 
itative behavior. This is particularly pronounced in the high density regime as 
dealt with in the next section. 

In Fig. 3 we compare the theoretical curves for the pressure-tensor compo- 
nents Pxx and Pyy with computer simulations for a density p = 0.01 and find 
excellent agreement. As noted above, it is evident that the square- root singular- 
ity of P originates from Pyy. 

It is interesting to note that in the narrow- channel regime, Ly < 2, the velocity 
distribution is still isotropic and is given by the Maxwell-Boltzmann distribution 
as is shown in Fig. 4. This supports the use of Eq. (11) for the average momentum 
transfer at a collision even for Ly <2. 

The singular behavior of the pressure is due to a sharp increase of the collision 
frequency as one approaches Ly = 2 from below. As can be seen in Fig. 2 the 
maximum Lyapunov exponent is reminiscent of this singularity as it is closely 
related to the collision frequency. At this low density this singularity is not clearly 
pronounced. However as will be shown below, at high densities the existence of 
a cusp at Lj^ = 2 is evident. The increase of the collision frequency is due 
to arresting configurations in which a pair of disks become trapped due to the 
boundaries. 



6 



> 



BB iiii H fliiflfliili fl ii 




0.5 



0.4 



0.3 



- 0.2 



0.1 



FIG. 2. Theoretical and simulation results for iV = 20 particles and a density 
p = 0.01 for the periodic-boundary case. Shown is the dependence of the maximum 
Lyapunov exponent Ai (right vertical scale), and of the equation of state {Pv/T) — 1 
(left vertical scale, full squares), and the rescaled single-particle collision frequency 
1^2 = (\/7r/4)t'2 (left vertical scale, open squares) on the channel width Ly. The full 
line is computed from Eq. (7). 
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FIG. 3. Comparison of the L^-dependence of the theoretical curves for {Pv/T) — 1, 
(Pxxv/T) — 1, and (Pyyv/T) — 1 (hues) with numerical simulations (points) for a system 
with 20 disks at a density p = 0.01. 
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III. HIGH DENSITY REGIME 

For hard-disk systems at high density one has to rely only on numerical simula- 
tions. In Fig. 5 we present the pressure, the collision frqucncy and the maximum 
Lyapunov exponent as a function of Ly at a density p = 0.4. It is observed that 
the singularity taking place at = 2 persists at high density in all three curves. 
According to Eq. (12) the curves for the collision freuqency and for Pv/T — 1 are 
proportional with a propertionailty constant independent of the density. This is 
evidently confirmed in Fig. 5. For the special case of only N = 2 particles in 
a square periodic box, the singularity at Ly — 2 appears at the critical density 
Pc = 0.5 and has recently been observed [7]. 

When the density is increased beyond 0.5, a broad peak emerges at some 
Ly < 2m the curves for the pressure and the maximum Lyapunov exponent. This 
is demonstrated in Fig. 6. This feature is attributed to the pronounced short- 
range crystalline-like order which is induced by the narrow-channel boundary 
conditions. Note that this order is not the natural triangular lattice of the system 
in the bulk. A typical microscopic configuration corresponding to this structure 
for a density 0.9 and Ly = 1.5 is shown in Fig. 7. The particles are arrested 
and cannot travel across the system, which has been referred to as the localized 
regime [7]. 

We have computed also the two pressure-tensor components P^x and Pyy for 
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FIG. 5. Simulation results for N = 20 particles and a density N/V = 0.4 for the 
periodic-boundary case. Shown is the dependence of the maximum Lyapunov expo- 
nent, Ai (right vertical scale), of the equation-of-state function (Pv/kT) — 1 (full cir- 
cles with label P, left vertical scale), and of the rcscalcd singlc-particlc collision fre- 
quency, u'2 = [(7r/A:r)^/^/4] U2 (crosses, left vertical scale) on the channel width Ly. 
v'2 agrees well with P. The dashed curves labelled Pxx and Pyy give the contributions 
(Pxxv/kT) — 1 and {Pyyv/kT) — 1 of the respective pressure-tensor components. 




a system with a density p = 0.8. As can be seen in Fig. 8, the two components 
are significantly different both above and below Ly = 2. The singularity of the 
pressure at Ly = 2 is evidently caused by Pyy, as is also the case for low densities. 

All the examples so far are for the case of periodic boundaries in both the x 
and y directions. A qualitatively similar result is obtained, if the long sides of the 
channel are elastically refiective, whereas the short sides remain periodic. This is 
demonstrated in Fig. 9 for = 100 particles at a density of N/V = 0.4, where 
P, Xi and u are shown as a function of the box width Ly. The maxima are not 
as sharp as in Fig. 5 for the periodic case, but are as pronounced. 

In order to get some insight into the dynamics of such systems it is usually 
useful to study the diffusion coefficient. It is well known that this coefficient does 
not exist in two-dimensional systems in the thermodynamic limit [15]. However, 
it is possible to compute the mean-squared displacement for a finite system and to 
extract effective diffusion constants and Dy in x and y directions, respectively. 
They are obtained from fits to the linear growth of (Aa;^) and (Ay^) for times 
t > 40. The results for a density p = 0.8 are shown in Figure 10. As expected 
D^. vanishes for Ly < 2, while Dy is non-zero. In fact, Dy increases with Ly 
approaching unity as a result of the periodic-boundary restrictions disfavoring 
momentum exchange in the y direction. Another intersting feature of Fig. 10 
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FIG. 7. A typical microscopic configuration for a density p = 0.9 and a cliannel widtli 
Ly = 1.5. It resembles the closest-dense packing consistent with the periodic bound- 
aries. It minimizes the mean free path and, hence, maximizes the collision frequency 
and the pressure. 
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FIG. 8. Simulation results for the pressure curves [Pv/T) — 1, {Pxxv/T) — 1, and 
{Pxxv/T) — 1 as a function of the channel width Ly. The system consists of A'^ = 20 
particles at a density p = 0.8. 
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FIG. 9. Simulation results for AT = 20 particles and a density N/V = 0.4 in a 
rectangular box with partly-reflecting boundaries (see the main text). Shown is the 
dependence of the maximum Lyapunov exponent, Ai (right vertical scale), of the equa- 
tion-of-state function {Pv/kT) — 1 (full circles with label P, left vertical scale), and 
of the rescaled single-particle collision frequency, 1^2 = 1^2 (crosses, left 

vertical scale) on the channel width Ly. v'^ agrees well with P. The dashed curves 
labelled P^x and Pyy give the contributions {PxxvjkT) — 1 and (Pyyv/kT) — 1 of the 
respective pressure-tensor components. 




FIG. 10. Effective diffusion constants (sec the main text) in x and y directions for 
disks in a narrow-channel with periodic boundaries at a density p = 0.8. 
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FIG. 11. Positive branches of Lyapunov spectra for 100-disk systems with periodic 
boundaries. The box width, Ly, is indicated by the labels. The density, p = 0.4, 
and the Volume, L^Ly, are fixed. The Lypunov exponents are defined only for integer 
values of the index I. In the inset the regime supporting Lyapunov modes is magnified. 
There, Ly increases from 1.2 (bottom) to 2.8 (top) in steps of 0.4. For details we refer 
to the main text. 

is the vanishing of and Dy for Ly > 2.7. This is due sohd-hke triangular 
configurations characteristic of the bulk at this high density. 

IV. SOME REMARKS ON LYAPUNOV SPECTRA 

In Fig. 11 we plot the Lyapunov spectra for a 100-disk system in a narrow 
periodic box. All spectra are for a density p = 0.4, which corresponds to a fluid. 
Since the volume, L^Ly, is fixed, an increase of the channel width Ly, varied in the 
figure in steps of 0.4, decreases L^ accordingly. Ly is indicated by the labels. Each 
spectrum consists of 4N Lyapunov exponents, of which only the positive branch, 
{A; > 0; / = 1, 2, . . . , 2A^}, is shown. The spectra are only defined for integer 
values of the index, /, which labels the exponents according to size. The lines are 
only drawn for clarity. According to the equilibrium formulation of the conjugate- 
pairing rule, the negative branch of the spectra, {A; >0;l — 2N + 1, . . . , 4A^}, is 
the mirror image of the positive branch [16]. 

The inset of Fig. 11 provides a magnified view of the small exponents. They 
turn out to be degenerate, which gives rise to a step-like appearance of the spec- 
tra [4]. The tangent-space perturbations associated with these exponents are 
collective, wave-like fields defined over the simulation box, and are referred to as 
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FIG. 12. The Lyapunov spectra plotted by the smooth lines are for 100-disk systems 
with periodic boundaries, where the. density, p = 0.9, and the volume, L^Ly, are 
fixed. The resulting box widths, Ly, are indicated by the labels. The dash-dotted 
curve labeled "solid" is for 100 disks in a periodic channel with an aspect ratio 
Ly/Lx = \/?>/{N/2) required for a triangular lattice. The density p is also given. 
For details we refer to the main text. Only the positive branches of the spectra are 
shown. 

Lyapunov modes [2,3,17]. The multiplicity is determined by the intrinsic sym- 
metries of the Hamiltonian and the boundary conditions, which give rise to the 
conserved quantities, energy and momentum. A complete classification in terms 
of transverse (T), longitudinal (L), and momentum (P) modes is given in Ref. 
[20]. Theoretical attempts have been made to interpret the modes in terms of 
fluctuating hydrodynamics [21-25]. 

We observe that a particular exponent, say the smallest positive, increases 
with Ly and, hence, with the wave number, k = 2t{/Lx = 2iTLy/V, belonging 
to this mode. This is referred to as a "dispersion relation" [2,3]. Thus, A2Af-3 
increases monotonously with Ly. However, as Fig. 11 shows, this proportion- 
ality does not hold for the large exponents. There is a crossover of the spectra 
reversing the sequence of A near Ly = 2. This is mainly a consequence of the 
collision frequency. The perturbations for the large exponents are found to be 
localized in space with only a small fraction of the particles contributing to the 
large exponents at any instant of time. The active zone moves around in space, 
restoring homogeneity on average. 

The smooth lines in Fig. 12 show Lyapimov spectra for analogous 100-disk 
systems for a density p = 0.9, where the channel width is indicated by the la- 
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bels. The case Ly = 1.5 corresponds to the arrested structure of Fig. 7 and is 
of particular interest. In Fig. 12 we compare it with the dash-dotted spectrum 
of 100 disks in a periodic channel with an aspect ratio Ly/Lx = V3/{N/2). For 
Ly = 1.8258 and N = 100, the density is p = 1.039, and the particles form 
a triangular lattice characteristic of a sohd in the bulk. Although the nearest 
neighbor separations and, hence, the collision frequencies, are nearly the same 
for both systems, the Lyapunov spectra are distinctly different. The maximum 
exponent, Ai, for the arrested structure is significantly larger than for the trian- 
gular system, whereas the behavior of the smaller exponents is just the reverse. 
This example demonstrates that the Lyapunov spectra are sensitive to structural 
details and are potentially a useful tool. But it is fair to say that no theory exists 
at present to interpret such spectra in full detail. 

V. DISCUSSION 

In this work the thermodynamic and dynamical behaviour of hard disks in a 
narrow channel is analysed using analytical and numerical approaches. The main 
results are for periodic boundary conditions in both x and y directions. Related 
studies of reflecting boundaries parallel to the channel axis reveal that similar 
features exist there as well. It is found that the pressure, collision frequency, 
maximum Lyapunov exponent, and Kolmogorov- Sinai entropy curves exhibit a 
singularity for a channel width equal to twice the disk diameter. For low densities 
this singularity for P and 1/2 is well understood by a virial expansion as it already 
shows up in the second virial coefficient. The singularity is not the result of a 
collective behavior and, in fact, is present for systems with a finite number of 
particles. Thus, it is not a genuine phase transition. It is rather a result of a 
rapid change in the available phase space taking place close to = 2. Similar 
but less pronounced singularities are expected to show up at = 3 and larger 
integers at higher virial coefficients. 

The singularity for the maximum exponent and for the KS-entropy is most 
pronounced for higher densities. This is a consequence of the fact that these quan- 
tities are closely related to the coUision frequency [11]. It has been established 
that in bulk systems these quantities exhibit a maximum at a phase transition 
[18,19]. The maximum Lyapunov exponent is a measure of the fastest dynamical 
events taking place in a system which are localized processes in space [3]. Thus, 
its behavior near the transition point is not related to a diverging length scale at 
the transition. Whereas in an ordinary phase transition in the bulk the enhanced 
collision frequency - and, thus, of Ai - is due to the emergence of a new structure, 
here it is a consequence of the constraints imposed by the boundaries. 

We expect the singularity found for Ly = 2 to be present also for soft disks. 
In this case, however, the singularity is likely to be weaker than the square-root 
singularity found in the pressure curve of hard disks. Similar features found in 
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the present study should be present in three-dimensional channels as well. It 
would be of interest to study these problems in more detail. 
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